{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modify the region 12 outlines in the RGI Regions files "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a cluster of glaciers South of the current extent of the region and subregion polygons. There are no regions below them, so there shouldn't be much of an issue in simply updating the geometry a little bit.\n",
    "\n",
    "Furthermore, the data type of the RGI_CODE attribute in the region file is int. For conistency with the RGI files, it should be str. We change this as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_dir = 'RGI62'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from oggm import utils\n",
    "import shapely.geometry as shpg\n",
    "import geopandas as gpd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/mowglie/Documents/git/rgitools/notebooks/rgi_60_corrections/RGI62/00_rgi62_regions'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out_dir = os.path.abspath(os.path.join(out_dir, '00_rgi62_regions'))\n",
    "utils.mkdir(out_dir, reset=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(32.0, 31.0, 54.0, 45.0)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rgi_reg = gpd.read_file(os.path.join(utils.get_rgi_dir(version='60'), '00_rgi60_regions', '00_rgi60_O1Regions.shp'))\n",
    "poly = rgi_reg.loc[rgi_reg.RGI_CODE == 12].iloc[0].geometry\n",
    "poly.bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's go down to 30° South instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = poly.exterior.xy\n",
    "ny = np.where(np.isclose(y, 31), 30, y)\n",
    "new_poly = shpg.Polygon(np.array((x, ny)).T)\n",
    "rgi_reg.loc[rgi_reg.RGI_CODE == 12, 'geometry'] = new_poly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Glacier 1\n",
    "# -129.7243, 54.4568 -> -129.7, 54.4492\n",
    "\n",
    "# Glacier 2\n",
    "# -126.9421, 57.7575 -> -126.9752, 57.7559\n",
    "# -126.9589, 57.7462 -> -126.9590, 57.7505"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in (1, 2):\n",
    "    poly = rgi_reg.loc[rgi_reg.RGI_CODE == i].iloc[0].geometry\n",
    "    x, y = poly.exterior.xy\n",
    "    x, y = np.array(x), np.array(y)\n",
    "    \n",
    "    ox, oy = -129.7243, 54.4568\n",
    "    nx, ny = -129.7, 54.4492\n",
    "    pok = np.argmin((x - ox)**2 + (y - oy)**2)\n",
    "    x[pok] = nx\n",
    "    y[pok] = ny\n",
    "    \n",
    "    ox, oy = -126.9421, 57.7575\n",
    "    nx, ny = -126.9752, 57.7559\n",
    "    pok = np.argmin((x - ox)**2 + (y - oy)**2)\n",
    "    x[pok] = nx\n",
    "    y[pok] = ny\n",
    "    \n",
    "    ox, oy = -126.9589, 57.7462\n",
    "    nx, ny = -126.9590, 57.7505\n",
    "    pok = np.argmin((x - ox)**2 + (y - oy)**2)\n",
    "    x[pok] = nx\n",
    "    y[pok] = ny\n",
    "    \n",
    "    rgi_reg.loc[rgi_reg.RGI_CODE == i, 'geometry'] = shpg.Polygon(np.array((x, y)).T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change type\n",
    "rgi_reg['RGI_CODE'] = [str(s) for s in rgi_reg.RGI_CODE]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg.to_file(os.path.join(out_dir, '00_rgi62_O1Regions.shp'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check\n",
    "rgi_reg = gpd.read_file(os.path.join(out_dir, '00_rgi62_O1Regions.shp'))\n",
    "assert rgi_reg.RGI_CODE.dtype == 'O'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(34.0, 32.0, 53.0, 42.0)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rgi_reg = gpd.read_file(os.path.join(utils.get_rgi_dir(version='60'), '00_rgi60_regions', '00_rgi60_O2Regions.shp'))\n",
    "poly = rgi_reg.loc[rgi_reg.RGI_CODE == '12-02'].iloc[0].geometry\n",
    "poly.bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = poly.exterior.xy\n",
    "ny = np.where(np.isclose(y, 32), 30, y)\n",
    "new_poly = shpg.Polygon(np.array((x, ny)).T)\n",
    "rgi_reg.loc[rgi_reg.RGI_CODE == '12-02', 'geometry'] = new_poly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in ('01-06', '02-02'):\n",
    "    poly = rgi_reg.loc[rgi_reg.RGI_CODE == i].iloc[0].geometry\n",
    "    x, y = poly.exterior.xy\n",
    "    x, y = np.array(x), np.array(y)\n",
    "    \n",
    "    ox, oy = -129.7243, 54.4568\n",
    "    nx, ny = -129.7, 54.4492\n",
    "    pok = np.argmin((x - ox)**2 + (y - oy)**2)\n",
    "    x[pok] = nx\n",
    "    y[pok] = ny\n",
    "        \n",
    "    rgi_reg.loc[rgi_reg.RGI_CODE == i, 'geometry'] = shpg.Polygon(np.array((x, y)).T)\n",
    "\n",
    "for i in ('01-06', '02-03'):\n",
    "    poly = rgi_reg.loc[rgi_reg.RGI_CODE == i].iloc[0].geometry\n",
    "    x, y = poly.exterior.xy\n",
    "    x, y = np.array(x), np.array(y)\n",
    "      \n",
    "    ox, oy = -126.9421, 57.7575\n",
    "    nx, ny = -126.9752, 57.7559\n",
    "    pok = np.argmin((x - ox)**2 + (y - oy)**2)\n",
    "    x[pok] = nx\n",
    "    y[pok] = ny\n",
    "    \n",
    "    ox, oy = -126.9589, 57.7462\n",
    "    nx, ny = -126.9590, 57.7505\n",
    "    pok = np.argmin((x - ox)**2 + (y - oy)**2)\n",
    "    x[pok] = nx\n",
    "    y[pok] = ny\n",
    "    \n",
    "    rgi_reg.loc[rgi_reg.RGI_CODE == i, 'geometry'] = shpg.Polygon(np.array((x, y)).T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_reg.to_file(os.path.join(out_dir, '00_rgi62_O2Regions.shp'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
